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Mechano-chemical spinodal decomposition: A phenomenological theory of phase 
transformations in multi-component, crystalline solids 

Shiva Rudraraju , Anton Van der Ven & Krishna Garikipati 

We present a phenomenological treatment of diffusion-driven martensitic phase transformations 
in multi-component crystalline solids that arise from non-convex free energies in mechanical and 
chemical variables. The treatment describes diffusional phase transformations that are accompa¬ 
nied by symmetry breaking structural changes of the crystal unit cell and reveals the importance 
of a mechano-chemical spinodal, defined as the region in strain-composition space where the free 
energy density function is non-convex. The approach is relevant to phase transformations wherein 
the structural order parameters can be expressed as linear combinations of strains relative to a high- 
symmetry reference crystal. The governing equations describing mechano-chemical spinodal decom¬ 
position are variationally derived from a free energy density function that accounts for interfacial 
energy via gradients of the rapidly varying strain and composition fields. A robust computational 
framework for treating the coupled, higher order diffusion and nonlinear strain gradient elasticity 
problems is presented. Because the local strains in an inhomogeneous, transforming microstruc¬ 
ture can be finite, the elasticity problem must account for geometric nonlinearity. An evaluation of 
available experimental phase diagrams and first-principles free energies suggests mechano-chemical 
spinodal decomposition should occur in metal hydrides such as ZrH2-2c- The rich physics that 
ensues is explored in several numerical examples in two and three dimensions and the relevance of 
the mechanism is discussed in the context of important electrode materials for Li-ion batteries and 
high temperature ceramics. 
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I. INTRODUCTION AND FORMULATION 

Spinodal decomposition is a continuous phase trans¬ 
formation mechanism occuring throughout a solid that 
is far enough from equilibrium for its free energy den¬ 
sity to lose convexity with respect to an internal degree 
of freedom. The latter could include the local composi¬ 
tion as in classical spinodal decomposition described by 
Cahn and Hilliard [ 1 ], or a suitable non-conserved order 
parameter as in Allen and Cahn’s theory for spinodal 
ordering [ 2 ]. A key requirement for continuous transfor¬ 
mations is that order parameters can be formulated to 
uniquely describe continuous paths connecting the var¬ 
ious phases of the transformation. These phases then 
correspond to local minima on a single, continuous free 
energy density surface in that order parameter space. For 
classical spinodal decomposition inside a miscibility gap, 
all phases have the same crystal structure and symmetry, 
and the order parameter is simply the local composition. 
The existence of a single, continuous free energy density 
surface for all phases participating in a transformation 
implies, by geometric necessity, the presence of domains 
in order-parameter space where the free energy density is 
non-convex. Reaching those domains through supersat¬ 
uration (by externally varying temperature or composi¬ 
tion) makes the solid susceptible to a generalized spinodal 
decomposition. 

Many important multi-component solids undergo 
phase transformations that couple diffusional redistribu¬ 
tion of their components with a structural change of the 
crystallographic unit cell. One prominent example is the 
decomposition that occurs when cubic yttria-stabilized 


zirconia Zri_cYc02-c/2 is quenched into a two-phase 
equilibrium region between tetragonal Zri_cYc02-c/2 
having low Y composition and cubic Zri_cYc 02 -c /2 hav¬ 
ing a high Y composition. Another occurs in Li battery 
electrodes made of spinel LicMn204. Discharging to low 
voltages causes the compound to transform from cubic 
LiMn204 to tetragonal Li2Mn2 04 through a two-phase 
diffusional phase transformation mechanism. As with 
simple diffusional phase transformations, these coupled 
diffusional/martensitic phase transformations can either 
occur through a nucleation and growth mechanism, or, if 
certain symmetry requirements are met, through a con¬ 
tinuous mechanism due to an onset of an instability with 
respect to composition and/or a structural order param¬ 
eter. 

Here we present a treatment of coupled diffu¬ 
sional/martensitic phase transformations triggered by 
instabilities with respect to both strain and composi¬ 
tion. These phase transformations are characterized by 
a mechano-chemical spinodal that is defined as a non- 
convex region of the free energy density function in the 
strain-composition space. The possibility of a mechano- 
chemical spinodal decomposition is motivated by recent 
first-principles studies of martensitic transformations in 
transition metal hydrides where a high temperature cu¬ 
bic phase is predicted to display negative curvatures with 
respect to strain, thus making strain a natural order pa¬ 
rameter to distinguish the cubic parent phase from its 
low temperature tetragonal daughter phases [ 33 ]. In ad¬ 
dition, the coupling with composition degrees of freedom 
(e.g. through the introduction of hydrogen vacancies in 
the metal hydrides) allows for the possibility that the free 
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energy also exhibits a negative curvature with respect to 
the composition. Mechano-chemical spinodal decomposi¬ 
tion, therefore, is a phenomenon that is likely present in 
many multi-component materials but has to date been 
overlooked as a mechanism by which a high symmetry 
phase can decompose martensitically and through dif- 
fusional redistribution upon quenching into a two phase 
regime. Structural phase transformations in solids driven 
by instability with respect to an internal shuffle of the 
atoms within the unit cell have been treated rigorously 
in the literature with coupled Cahn-Hilliard and Allen- 
Cahn approaches [ 5 , 6, 7 , 8]; but mechano-chemical spin¬ 
odal decomposition is fundamentally different and neces¬ 
sitates a coupled treatment of both the strain and com¬ 
position instabilities. 

Our treatment is based on a generalized, Landau- 
type free energy density function that couples strain 
and composition instabilities. The governing equations 
of mechano-chemical spinodal decomposition, obtained 
by variational principles, generalize the classical equa¬ 
tions of the Cahn-Hilliard formulation [1], and of nonlin¬ 
ear gradient elasticity [ 9 , 11 ] by coupling these systems. 
The ability to solve this complex, nonlinear, strain- and 
composition-gradient-driven, mechano-chemical system 
for sufficiently general initial and boundary value prob¬ 
lems in two and three dimensions also has been lacking 
heretofore. We introduce the computational framework 
to obtain such solutions in general, three-dimensional 
solids. This newfound capability allows us to then rein¬ 
force our discussion of the phenomenology with dynamics 
predicted by the numerical solutions. 


A. The mechano-chemical spinodal in two 
dimensions 

For accessibility of the arguments, we first consider 
the two-dimensional analogue of the cubic to tetrago¬ 
nal transformation: the square to rectangle transforma¬ 
tion. The high-symmetry square lattice will serve as the 
reference crystal relative to which strains are measured. 
Lower symmetry lattices that can be derived from the 
square lattice by homogeneous strain include the rect¬ 
angle, the diamond and lattices where there are no con¬ 
straints on the cell lengths and their angles. 

We use the composition c, varying between 0 and 
1, as our order parameter for the chemistry of our bi¬ 
nary two-dimensional solid. Symmetry breaking struc¬ 
tural changes are naturally described by the strain rel¬ 
ative to the high symmetry square lattice. The elastic 
free energy density is also a function of strain. In both 
contexts, strain will in general be of finite magnitude. 
Following Barsch & Krumhansl [ 9 ], we use the Green- 
Lagrange strain tensor, relative to the square lattice. 
Rotations are exactly neutralized in this strain measure; 
for any rigid body motion, E = 0 (Supporting Informa¬ 
tion). In two dimensions, the relevant strain components 
are Eii^E22 and E12 = E21. However, it is more conve- 



FIG. 1: (a) Square to rectangle transition (2D) in 

reparametrized strain space, (b) Cubic to tetragonal transi¬ 
tion (3D) in reparametrized strain space (Equation[1]). Lat¬ 
tice vectors are labelled by the corresponding coordinate 
directions xi,x 2 ,x 3 and the corresponding lower symmetry 
phases are labelled 1, 2, 3. The deformations shown in the 
sub-figures are area/volume preserving, respectively. 


nient to use linear combinations of these components that 
transform according to the irreducible representations of 
the point group of the high symmetry square lattice. In 
two dimensions, these include ei = {En + F^22)/v^, 
62 = {Ell — E22)IV^ and cq = \/ 2 Ei 2 . Here, ei and 
66 reduce to the dilatation and shear strain, respectively 
in the infinitesimal strain limit. The strain measure 62 
uniquely maps the square lattice into the two rectangu¬ 
lar variants [Figure la]: positive and negative 62 generate 
the rectangles elongated along the global Xi and X2 di¬ 
rections, respectively. The equivalence of the rectangular 
variants under the point group symmetry of the square 
lattice (62 = 0) restricts the free energy density to even 
functions of 62. 

If the crystalline solid has multiple chemical species, 
its free energy density dependence on 61, 62 and 66 can 
change with composition, c. Figure 2 a illustrates a free 
energy density surface, ^(6,62), for a binary solid that, 
at higher temperature, forms a solid solution having 
square symmetry. In this case ^ is convex, a condi¬ 
tion made precise by specifying positive eigenvalues of 
its Hessian matrix over the {6,62} space. Additionally, 
^(6, 0) is a minimum with respect to 62 for fixed 6, mak¬ 
ing the square phase stable for all c at this tempera¬ 
ture. However, at a lower temperature, ^ may lose con¬ 
vexity, inducing the notion of a mechano-chemical spin¬ 
odal region. We define it as the domain in {6,62} space 
over which the Hessian matrix admits negative eigen¬ 
values, as illustrated in Figure 2 b. Here, we focus on 
the conditions jdc? < 0 and jde^ < 0 . The 
square phase (62 = 0) remains stable at high compo¬ 
sition ( 6 ^ 1 ) with positive eigenvalues of the Hessian 
(Figure 2 b). But it is mechanically unstable at low com¬ 
position, with ldc\ < 0 for (6,62) ^ (0,0), and has 
two symmetrically equivalent, stable, rectangular phases 
with ldc\ > 0 at (e, 62) = ( 0 , Esq). While not shown 
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in Figure 2 b we assume that ^ is convex with respect to 
ei and Figure 2 c illustrates a schematic temperature- 
composition phase diagram consistent with the free en¬ 
ergy densities of Figure 2 a-b. The square/cubic phase 
forms at high composition or at high temperature; the 
rectangle/tetragonal phase a, forms at low composition 
and low temperature. A large two-phase region separates 
them. 

Zuzek et al. [ 32 ] have obtained phase diagrams ex¬ 
perimentally that are topologically equivalent to Figure 
2 c for the ZrH2-2c system. Recent first principles cal¬ 
culations [ 33 ] have demonstrated the existence of a me¬ 
chanical instability that exists in this system at low c via 
non-convexity with respect to strain. Furthermore, the 
two-phase coexistence at low temperature also implies 
non-convexity with respect to composition as we have 
demonstrated in supporting information. Figures 2 a and 
b therefore represent a two-dimensional analogue of the 
free energy for such systems. 




FIG. 2: (a) Free energy density for the 2D formulation at 
high temperature, and (b) at low temperature, showing the 
mechano-chemical spinodal along with two energy minimizing 
paths (brown lines) and their corresponding minimum energy 
strained structures, (c) Temperature-composition phase dia¬ 
gram. 


or negative 62, driven towards a local minimum at con¬ 
stant c. These deformations due to the mechanical insta¬ 
bility will happen like many martensitic transformations 
where a mix of symmetrically equivalent rectangular vari¬ 
ants coexist to minimize macroscopic strain energy. For 
finite but moderate strain, the transformation could pro¬ 
ceed coherently, even if the two symmetrically equivalent 
rectangular variants coexist. We neglect non-essential 
complexities of this process and assume that, instantly 
upon quenching, finitely sized neighborhoods of the solid 
deform homogeneously into one of these rectangular vari¬ 
ants at the original composition. 

The solid also becomes susceptible to uphill diffusion 
because jdc? < 0 implies a negative diffusion coef¬ 
ficient. However, it does not occur at constant 62, since 
the valleys traversing the local minima, Ide2 = 0, 
between the square lattice at c = 1 and the rectangular 
lattices at c = 0 span intervals of negative and positive 
62- Mechano-chemical spinodal decomposition sets in. 
Composition modulations are amplified: high c regions 
strive to be more square (point C) while low c regions 
strive to be more rectangular (points D or D'). How¬ 
ever, since coherency is maintained, some neighborhoods 
in the solid will be frustrated from attaining strains that 
ensure minima, 9^/9e2 = 0, for local values of c. Co¬ 
herency strain-induced free energy penalties arise to alter 
the driving forces for purely chemical spinodal decompo¬ 
sition. Movie S 5 in the supporting information shows the 
evolution of the state (c, 62) of the material points on the 
free energy manifold ^. 


B. Mathematical formulation: three dimensions 

Armed with the insight conveyed by the two- 
dimensional study, we next lay out the three-dimensional 
treatment, in which setting the considered phase trans¬ 
formations proceed via lattice deformation and diffusion. 
The arguments are made more concrete by considering 
the general mathematical form of the free energy den¬ 
sity. 


1. Strain order parameters 


Consider our model binary solid, annealed at high tem¬ 
perature Th, to form a solid solution in the square phase, 
13 . Its state is at point A in Figure 2 a, with 62 = 0 . It 
is then quenched into the two-phase region (Figure 2 c) 
with free energy density at point B in Figure 2 b. For 
a quench at sufficiently high rate, the dimensions of the 
square lattice, controlled by strains ei, 62 and ee, and 
the composition remain momentarily unchanged. How¬ 
ever, since the state at point B satisfies jde^ < 0 
and jdc? < 0, there exist thermodynamic driving 
forces for segregation by strain and composition within 
the mechano-chemical spinodal. 

Diffusion being substantially slower than elastic relax¬ 
ation, the solid immediately deforms to either positive 


The strain measures ei-ee are first redefined using the 
full Green-Lagrange strain tensor components in three 
dimensions: 


ei — + E22 + ^33)5 ^2 — — ^22)5 

^3 = + ^22 — 2F^33) 


65 


= 


31 , 


64 = V2E23 = V^^ 32 , 

Cq = ^/2Ei2 = ^/2E2\ ( 1 ) 


In the limit of infinitesimal strains, ei describes the di¬ 
latation, while 64, 65 and ee reduce to shears. The point 
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group operations of the cubic crystal leave ei invariant 
since it is the trace of while collecting each of the sub¬ 
sets { 62 , 63 } and { 64 , 65 , 66 } into a symmetry-invariant 
subspace whose elements transform into each other. The 
measures 62 and 63 are especially suited as order parame¬ 
ters to describe cubic to tetragonal distortions. All three 
tetragonal variants that emerge from the cubic reference 
crystal can be uniquely represented by these two mea¬ 
sures. See Figure lb where tetragonal distortions of the 
cubic crystal along the Xi, X 2 and X 3 axes have been 
labelled as 1 , 2 and 3, respectively. Deviations from the 
dashed lines in the 62-63 space of Figure lb correspond 
to orthorhombic distortions of the cubic reference crys¬ 
tal. This role of 62 and 63 as structural order parameters 
to denote the degree of tetragonality, and to distinguish 
between the three tetragonal variants, complements their 
fundamental purpose as arguments of the elastic free en¬ 
ergy density. 



FIG. 3: Mechano-chemical spinodal for the 3D formulation 
depicted by contours of the free energy manifold along the 62 - 
63 and es-c planes. The three energy minimizing paths (brown 
lines) and their corresponding energy minimizing strained 
structures and are also shown. 


2. The mechano-chemical spinodal in three dimensions 

For brevity we write the strains as e = { 61 ,..., ce}. 
We introduce further phenomenology by specifying that 
^{c,e) is only one part of the total free energy density. 
It is a homogeneous contribution whose composition 
and strain dependence cannot be additively separated. 
Figure 3, for example, illustrates contour plots of ^(e, e) 
on the 62-63 and 63-6 planes for a binary solid having 
a temperature-composition phase diagram similar to 
that of Figure 2c. To generate the tetragonal a phase 
as illustrated in that diagram, the homogeneous free 
energy density as a function of strain must qualitatively 
follow the contours of the 62-63 plane at 6 = 0 in Figure 
3. Here, the tetragonal variants are local minimizers of 
^ in 62-63 space, with equal minima. In turn, to obtain 


the cubic /3 phase, ^(e, e) must follow the contours of 
the 62-63 plane at 6 = 1. Thus, ^ changes smoothly 
from convex with respect to 62 and 63 on the 6=1 plane 
to non-convex at 6 = 0 to form the three variants of the 
tetragonal a phase. Importantly, the planes at 6 = 0 
and 6=1 themselves must be minimum energy surfaces 
to represent the tetragonal a and cubic f3 phases, 
respectively. Movie S 8 in the supporting information 
shows the evolution of the state (e, 62 , 63 ) of the material 
points on the free energy manifold 

Gradient regularization of the free energy density: Fol¬ 
lowing van der Waals [10], and Cahn & Hilliard [1] in 
their treatment of non-uniform composition fields, we can 
extend the total free energy density beyond ^ by writ¬ 
ing it as a Taylor series, retaining terms that depend on 
the composition gradient, Ve. We extend this gradient 
dependence to the strain measures Ve, as did Barsch & 
Krumhansl [9] following Toupin [ 11 ]. (Also see Karatha 
and co-workers [12, 13], for treatments using infinitesi¬ 
mal and finite strain, respectively.) These gradient de¬ 
pendences appear in a non-uniform free energy density 
^(6, e, V6, Ve). Frame invariance of ^ and ^ is guar¬ 
anteed since the members of e are linear combinations of 
the tensor components of E. They also must be invari¬ 
ant under point group operations of the cubic reference 
crystal. 


3. Free energy functional 

The crystal occupies a reference (undeformed) config¬ 
uration D with boundary F. The total free energy, H, is 
the integral of the free energy density ^ ^ over the 

solid with boundary contributions included. Thus, H is 
a functional of the composition c and the displacement 
vector field n, from which the strains are derived (see 
Supporting Information): 



UiTidS. 


( 2 ) 


where traction vector component Ti is specified on the 
boundary subset F^. C F. Following the above authors 
we include only quadratic terms in the gradients, but in 
a generalization we also allow cross terms between Ve 
and Ve in the non-uniform contribution which can 
therefore be written as 


^(6, e, Ve, Ve) =^V6 • ^(e, e)V6 

a,3 

+ ^Vc-6»“(c,e)Ve„. (3) 

a 
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Here k is a symmetric tensor of composition gradient en¬ 
ergy coefficients, each 7*^^ = 1,...,6) is a tensor 

of strain gradient energy coefficients, and each 6 ^ is a 
tensor of the mixed, composition-strain gradient energy 
coefficients. Note that, in general, these coefficients will 
be functions of the local composition and strain. The 
point group symmetry of the cubic reference crystal im¬ 
poses constraints on the tensor components of k, 7^^ and 
9 ^ as well as on the form of 

While the gradient energies bestow greater accuracy 
upon the free energy description of solids with non- 
uniform composition and strain fields, they are essential 
at a more fundamental level if the homogeneous free en¬ 
ergy density is non-convex. At compositions that render 
^ non-convex, the absence of a gradient energy term will 
allow spinodal decomposition characterized by composi¬ 
tion fluctuations of arbitrary fineness, thus leading to 
non-unique microstructures—a fundamentally unphysi¬ 
cal result [ 14 ]. With k / 0, the composition gradient en¬ 
ergy penalizes the interfaces wherein composition varies 
rapidly between high and low limits. This ensures phys¬ 
ically realistic results, manifesting in unique microstruc¬ 
tures with a mathematically well-posed formulation. 

An essentially analogous situation exists with respect 
to the negative curvatures of in the 62-63 plane at 
low 6, which drive the cubic lattice to distort into the 
tetragonal variants corresponding to the three free en¬ 
ergy wells. Consider a solid with a homogeneous free 
energy density as in Figure 3 and a strain state lying be¬ 
tween the valleys in the 62-63 plane at 6 = 0. Absent the 
strain gradient energy, mechano-chemical spinodal de¬ 
composition would allow tetragonal variants of arbitrary 
fineness—an unphysical result, reflecting further math¬ 
ematical ill-posedness. Retention of the strain gradient 
energy, (7*^^ ^ 0) penalizes interfaces of sharply vary¬ 
ing strain between tetragonal variants to ensure phys¬ 
ically realistic results and unique microstructures from 
a mathematically well-posed formulation. This is well- 
understood in the literature that studies the formation 
of martensitic microstructures from non-convex free en¬ 
ergy density functions [ 16 - 18 ]. 


4^ Governing equations of non-equilibrium ehemistry 

The free energy for non-homogeneous composition and 
strain fields, Eq. [ 2 ], must be a minimum at equilibrium. 
The state of a solid out of equilibrium will evolve to re¬ 
duce the free energy n[6, n]. In formulating a kinetic 
equation for the redistribution of atomic species through 
diffusion, we are guided by variational extremization of 
the free energy to identify the chemical potential, fi. De¬ 
tails of this calculation appear as supporting information. 
The result follows: 


. _ dK_ 

ti =— -V • (kVc) + Vc • ^V6 

oc oc 



a,l3 


+ ^fVc-—Ve«-V-(6»“Ve„)j. ( 4 ) 


For solids where c tracks the composition of an intersti¬ 
tial element within a chemically inert host, such as Li 
in LicMn204, /i in Eq. [ 4 ] corresponds to the chemical 
potential of the interstitial element. If c tracks the com¬ 
position of a substitutional species, such as in alloys or 
on sublattices of complex compounds (e.g. the cation 
sublattice of yttria stabilized zirconia), /r is equal to the 
chemical potential difference between the substitutional 
species. 

The common phenomenological relation for the flux 
is J = —T(6, e)V//, where L is the Onsager transport 
tensor (see de Groot & Mazur [ 15 ]). For an interstitial 
species, L is related to a mobility [ 19 ], while it is a kinetic, 
intermixing coefficient for a binary substitutional solid 
[ 20 ]. Inserting the flux in a mass conservation equation 
yields the strong form of the governing partial differen¬ 
tial equation (PDF) for time-dependent mass transport. 
It is of fourth order in space due to the composition gra¬ 
dient dependence of /i in Equation [ 4 ]. See Supporting 
Information for strong and weak forms of this PDE. 


5. Governing equations of meehanieal equilibrium: Strain 
gradient elastieity 

Mechanical equilibrium is assumed since elastic wave 
propagation typically is a much faster process than dif- 
fusional relaxation in crystalline solids. Equilibrium is 
imposed by extremizing the free energy functional with 
respect to the displacement field. Standard variational 
techniques lead to the weak and strong forms of strain 
gradient elasticity. The treatment is technical, for which 
reason we restrict ourselves to the constitutive relations 
that are counterparts to the chemical potential equation 
[ 4 ] for chemistry. Coordinate notation is used for trans¬ 
parency of the tensor algebra, and summation is implied 
over repeated spatial index, /. Details appear in Sup¬ 
porting Information. The final form of the equations is 
complementary to Toupin’s [II], since our derivation is 
relative to the reference crystal, Gt. 

With the deformation gradient F being related to the 
Green-Lagrange strain as Ekl = \{FiKFiL — Srl)^ the 
first Piola-Kirchhoff stress tensor and the higher-order 
stress tensor respectively, are given by 
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PiJ 


BiJK 


E 


d{^ + dCa 

dca dFij 


^ dCaJ 

^ dcaj dFij 


E 


deaj 
dcaj dFij^K 


( 5 ) 

(6) 


The higher-order stress which is absent in classical, 
non-gradient elasticity [21] (and in earlier treatments of 
mechano-chemistry [ 4 , 22 ]) makes the strong form of gra¬ 
dient elasticity a fourth order, nonlinear PDE in space 
(Supporting Information). The first three-dimensional 
solutions to general boundary value problems of Toupin’s 
strain gradient elasticity theory at finite strains were re¬ 
cently presented by the authors [ 23 ]. 


II. RESULTS 



mesh deformation 


FIG. 4: Evolution of 2 D microstructure during outflux from 
the top and bottom surfaces of a solid under plane strain. 
Contours show strain 62 . Note the legend and correspond¬ 
ing square/rectangular variant crystal structures. The defor¬ 
mation and accompanying twinned microstructure are seen 
clearly in the distorted mesh. 


A. Two dimensional examples 

We first consider a two-dimensional solid to bet¬ 
ter visualize the microstructures that can emerge from 
mechano-chemical spinodal decomposition. Plane strain 
elasticity is assumed, for which ^13,^23,^33 = 0, giv¬ 
ing 64,65 = 0 , and 63 = 6 i/\/ 2 , reducing Equa¬ 
tions [ 1 - 6 ] to two dimensions. The discussion on two- 
dimensional mechano-chemical spinodal decomposition 
(Eigure 2) holds: as a particular parameterization of 
.^(6,61,62,66) we consider a regular solution model as 
a function of c at zero strain. At 6 = 0, .^(6,61,62,66) 
is double-welled in 62 corresponding to the two rectan¬ 
gular variants. The gradient energy is ^(Ve, V62) with 
a constant, isotropic k 7^ 0 , and constant 7^^ 7^ 0, all 
other gradient coefficients being zero. While the fullest 
possible complexity of the coupling is not revealed by 
these simplifications, the aim here is to present the essen¬ 
tial physics that is universal to mechano-chemical spin¬ 
odal decomposition, postponing details of the more com¬ 
plex couplings and tensorial forms to future communi¬ 
cations, where they will be derived for specific materials 
systems. See Supporting Information for specific forms 
of ^ (c, 61,62,66) and ^(V6, V62). 

Eigure 4 shows the evolution of microstructure over 
a 0.01 X 0.01 domain whose reference (initial) state has 
the square crystal structure and 6 = 1 . The displace¬ 
ment component Ui = 10“^ on the right boundary {Xi = 
0 . 01 ), with the remaining boundaries fixed (n = 0 ). An 
outward fiux is imposed on the top and bottom bound¬ 
aries causing a decrease in composition, 6, starting from 
the boundaries. As the composition falls, the homo¬ 
geneous free energy density ^ loses convexity and the 
state of the material (e, 62) enters the mechano-chemical 
spinodal along 62 = 0 (Eigure 2 b). The continuing out¬ 
ward flux first drives material near the top and bottom 
boundaries fully into the regime where the rectangular 
crystal structure is stable. As explained for Eigure 2b, 


the negative curvature jde^ < 0 creates thermody¬ 
namic driving forces that distort the square structure, 
shown in green, into rectangular variants, which form as 
red/blue laminae (all 62 values). The coexistence of the 
parent, square structure with the daughter, rectangular 
variants also can give rise to cross-hatched microstruc¬ 
tures that parallel the tweed microstructures described 
in the work of Karatha and co-workers [12, 13 ]. Movie 
S 9 in the supporting information shows the formation of 
such a microstructure. 

A laminar, micro-twinned microstructure develops as 
the two rectangular variants form, distinguished by the 
sign of 62 (see legend). Note that 62 remains continu¬ 
ous because of the penalization of strain gradients V62, 
although discontinuities can develop in V62, itself. Be¬ 
cause our numerical framework uses basis functions that 
are continuous up to their first derivatives (see Methods 
and Supporting Information), V62 is indeed discontinu¬ 
ous in the computations. The lamination accommodates 
the strain difference between the two rectangular variants 
to minimize the free energy. If strains were infinitesimal 
(| 6 iI, 1621, 1661 ^ 0 ) 62 would correspond to shear along 
directions that are rotated 7r/4 radians from the crystal 
axes. But the strains in these microstructures can be fi¬ 
nite (reflected in the range — 0.1 < 62 < 0.1 in Eigure 4 ) 
and involve large rotations to accommodate the rectangu¬ 
lar variants; this necessitates a finite-strain formulation, 
as shown in recent studies by Einel et al [ 41 ] compar¬ 
ing infinitesimal-strain and finite-strain formulations for 
polytwinned microstructure evolution in martensitic al¬ 
loys, and by Hildebrand & Miehe [8]. 

The micro-twinning and coherency strains are seen in 
the distorted mesh of discretization (inset of top cen¬ 
ter panel). The undeformed mesh cells are squares, and 
hence the numerical discretization strikingly delineates 
the kinematics of the cubic to rectangular transforma¬ 
tion, and highlights the rectangular twins as well as 
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the concentrated distortions along the twin boundaries. 
Movie S6 in the supporting information shows such mesh 
distortion and formation of the rectangular twins with a 
different form of the function ^. 

In some cases, the long-range nature of elasticity forces 
like rectangular variants to align even when separated 
by an untransformed square phase. This is first seen 
in the finger-like extensions of strain contours from the 
rectangular variants into the square phase in Figure 4 , 
followed by their alignment and eventual incorporation 
into laminae of the same variant (top right panel and its 
evolution shown as an inset). In this instance, although 
the micro-twins end at an invariant habit plane that is 
common with the untransformed square phase, the lattice 
parameters of the rectangular variants differ from those 
of the square phase, inducing elastic strain energy in the 
latter. The alignment lowers this strain energy, and also 
is seen in Figure 5 a. Note, however, that this is not a 
universal feature. Movie S2 shows instances where unlike 
variants come close to alignment. 



ll = l.O 11 = 0.1 11 = 0.01 


FIG. 5: Microstructure controlled by elastic gradient length 
scale parameter (le). Shown are the final contours of 62 for 
the simulation of (a) outflux from top and bottom surfaces, 
and (b) quenching. 

The fineness of laminae depends on the 

strain gradient elasticity length scale Iq ^ 

\j y^,/3=i,2.6(^^'^/^ea^e/3)^ as explored in 

Figure 5 . For Figure 5 a, the initial and boundary 
conditions are the same as in the example of Figure 
4 . Decreasing weakens the penalty on the strain 
gradient Ve2 across neighboring, unlike rectangular 
variants and allows more twin boundaries. Notably, 
self-similarity is not maintained between microstructures 
for different /g, even for the same initial and boundary 
value problem. We understand this to be the influence 
of elastic strain accommodation: To minimize the total 
free energy when the strain gradient penalty changes, 
the physics optimizes twin boundaries via laminations of 
different sizes as well as different patterns. Importantly 
however, crystal symmetry admits non-vanishing strain 


gradient energy coefficients beyond 7^^ 7^ 0 used in these 
simulations (Supporting Information). Furthermore, the 
composition dependence of ^ could be more complex 
than the simple regular solution model used here. Given 
the already strong effect of /g alone, we conjecture that 
varying these forms will have a significant influence on 
the resulting coherent twin microstructure. The proper 
form, while guided by crystal symmetry arguments 
must ultimately be determined experimentally or by 
first-principles statistical mechanical methods. 

The example in Figure 5 b further explores the influ¬ 
ence of /g. In this suite of computations, the unstrained 
material with an initially square microstructure, convex 
free energy density and composition having random fluc¬ 
tuations about c = 0.45 was quenched to a low tempera¬ 
ture and into the mechano-chemical spinodal. The local 
composition and strain evolve under the thermodynamic 
driving forces detailed in the context of Figure 2. We 
draw attention, once again, to the changing identity of 
rectangular variants, their shapes and sizes depending on 
/g. Also note the progressively finer lamination of rectan¬ 
gular domains with decreasing /g. Such studies suggest 
how dynamic mechano-chemical spinodal decomposition 
can lead to an atlas of microstructures, which in turn will 
determine material properties. Movies S 1 -S 4 in the sup¬ 
porting information show the evolution of some of these 
microstructures. 


B. A three dimensional study 

This final example displays the full three-dimensional 
complexity of microstructures resulting from the 
mechano-chemical spinodal. We persist with the above 
simplifications aimed at presenting the essential physics 
of mechano-chemical spinodal decomposition, and post¬ 
pone a full exploration of the coupling and anisotropies 
in coefficients to future work on specific materials sys¬ 
tems. The free energy density function used for the three- 
dimensional study appears as supporting information. 

Figure 6 shows the equilibrium microstructure that re¬ 
sults in a solid that is initially in the cubic phase, c = 1, 
subject to an outflux on all surfaces. The domain is a unit 
cube with displacement components U2,us = 0.01 on the 
boundary Xi = 1, with zero displacement, n = 0, on 
the boundary Xi = 0 . The cubic to tetragonal trans¬ 
formation takes place as c ^ 0 . Under these condi¬ 
tions all three tetragonal variants form, with an intri¬ 
cately interleaved microstructure for strain accommoda¬ 
tion. Note the finer microstructures and changing pat¬ 
tern for smaller /g. The inset shows the surface strain¬ 
ing around a corner, delineated by the distorted mesh 
lines. Observe the three, oriented, tetragonal variants 
formed by twinning deformation from the initial cubic 
structure. See Movie S 7 in the supporting information 
for a detailed view of the three-dimensional structure of 
these individual tetragonal variants. Movie SIO in sup¬ 
porting information shows other three-dimensional mi- 
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crostructures whose cross-sectional planes bear closer re¬ 
semblance to the plate-like structures predicted by two- 
dimensional computations. To the best of our knowledge, 
such computations of a cubic to tetragonal transforma¬ 
tion with twinned variants whose microstructure is con¬ 
trolled by nonlinear, strain gradient interface energies, 
have not been previously presented. 



FIG. 6: Microstructure observed in 3D for different values of 
the elastic gradient length scale parameter (le). The three 
tetragonal variants appear in blue (variant 1), yellow (variant 
2) and red (variant 3) for c = 0. The transformation strains 
are easily discerned in the distorted mesh. 


III. DISCUSSION 

Several kinetic mechanisms have been proposed to de¬ 
scribe the decomposition of a solid upon entering a two- 
phase region [ 42 , 43 ]. A commonly observed mechanism 
occurs through localized nucleation followed by diffu- 
sional growth that is mediated by the migration of in¬ 
terfaces separating the daughter phases from the parent 
phase. Nucleation and growth mechanisms are common 
when the phases participating in the transformation have 
little or no crystallographic relation to each other. The 
transformation then proceeds reconstructively, where the 
interphase interfaces are disordered or at best only semi- 
coherent. Numerical treatments of this mechanism rely 
on sharp-interface methods such as a level set approach 
[ 44 ]. 

Other kinetic mechanisms of decomposition involving 
some degree of coherency are also possible when the par¬ 
ticipating phases share sufficient crystallographic com¬ 
monality that order parameters can be defined describing 
continuous pathways connecting the parent and daugh¬ 
ter phases. The structural changes of the crystal that 
accompany these transformations can fall in one of two 
categories. One subclass of structural transformations is 
driven by an internal shuffie, where the atomic arrange¬ 
ment within the unit cell of the crystal undergoes a sym¬ 
metry breaking change. The unit cell vectors of the crytal 
may also undergo a change and lower the symmetry of the 


lattice, but this effect is secondary and in response to the 
atomic shuffie within the unit cell. The order parameters 
for the structural change therefore describe the atomic 
shuffles within the unit cell and are non conserved. The 
second subclass of structural transformations are driven 
by a symmetry reducing change of the lattice vectors of 
the crystal, with any atomic rearrangement of the ba¬ 
sis of the crystallographic unit cell occurring in response 
to the symmetry breaking changes of the unit cell vec¬ 
tors. The natural order parameters describing these tran¬ 
sitions are strains. The free energy in the first subclass 
will exhibit negative curvatures with respect to the non- 
conserved shuffie order parameters, while the free energy 
in the second subclass will have negative curvatures with 
respect to the relevant strain order parameters, as was 
recently predicted for the cubic to tetragonal transition 
of ZrH2 [ 33 ]. 

Decomposition reactions that require both diffusional 
redistribution and a structural change due to an internal 
shuffle have been treated successfully and rigorously with 
phase field approaches that combine Cahn-Hilliard with 
Allen-Cahn. The Cahn-Hilliard description accounts for 
the composition degrees of freedom while the Allen-Cahn 
component describes the evolution of the non-conserved 
shuffie order parameters required to distinguish the vari¬ 
ous phases of the transformation [ 5 - 8 ]. These treatments 
have included strain energy as a secondary effect serving 
only as a positive contribution to the overall free energy 
due to coherency strains. The approach is rigorous if the 
free energy density remains convex with respect to strain, 
with instabilities in the free energy appearing as a func¬ 
tion of concentration and the non-conserved shuffie order 
parameters. 

Here we have presented a mathematical formulation 
and an accompanying computational framework for de¬ 
composition reactions that combine diffusional redistri¬ 
bution with a structural change driven by a symmetry 
breaking strain of the crystallographic unit cell as op¬ 
posed to an internal shuffie within the unit cell. Strains 
play a primary role, explicitly serving as order parame¬ 
ters to distinguish variants of a daughter phase that has 
a symmetry subgroup-group relationship with its parent 
phase due to a structural change of the crystallographic 
unit cell. Crucial to the description is that the driving 
force for the formation of the daughter from a super¬ 
saturated parent phase emerges from an instability with 
respect to composition. The treatment therefore com¬ 
bines Cahn-Hilliard for composition with a description 
of martensitic transformations at fixed concentration in¬ 
troduced by Barsch and Krumhansl [ 9 ]. The existence 
of simultaneous instabilities with respect to strain and 
composition has not been considered as a mechanism to 
describe decomposition reactions upon quenching into a 
two-phase region. 

The possibility of such a mechano-chemical spinodal 
mechanism is motivated by recent first-principles studies 
of the cubic to tetragonal phase transformations of ZrH2, 
where the free energy of the high temperature cubic form 
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was predicted to become unstable with respect to strain 
upon cooling below the cubic to tetragonal second order 
transition temperature [ 33 ]. The ZrH2-2c hydride can ac¬ 
commodate large concentrations of hydrogen vacancies, 
c, and has a phase diagram that is topologically identical 
to that depicted in Figure 2 c, with a two-phase region 
separating a hydrogen-rich tetragonal form of ZrH2-2ca 
from a cubic form of ZrH2-2c/3 (with Co, < C/3). See Zuzek 
et ah [ 32 ].[ 45 ] To be consistent with the predicted free 
energies for stoichiometric ZrH2 (i.e. c=0) and the ex¬ 
perimental T versus c phase diagram with the form of 
Figure 2 c, the free energy of this hydride as a function 
of composition and strain (i.e. 62 and 63) should be sim¬ 
ilar to those depicted in Figures 2 a and b, and Figure 
3 . Decomposition upon quenching cubic ZrH2-2c into 
the two-phase region can therefore proceed through a 
mechano-chemical spinodal decomposition mechanism. 

Not only does the phenomenological description in¬ 
troduced in this work describe a new mechanism of 
decomposition that is qualitatively distinct from previ¬ 
ous combined Cahn-Hilliard and Allen-Cahn approaches, 
its numerical solution also proves to be substantially 
more complex due to contributions from strain gradi¬ 
ent terms. Indeed, even numerical solutions to general, 
three-dimensional boundary value problems of gradient 
elasticity at finite strain were not available until pre¬ 
sented recently by the authors [ 23 ]. Furthermore, as 
other authors have pointed out before, the use of strain 
metrics as order parameters makes a reliance on infinites¬ 
imal strains untenable due to the rigid rotations that ac¬ 
company the finite strains characterizing most structural 
transformations [8, 41 ]. The use of finite strain met¬ 
rics introduces geometric non-linearity into the problem, 
which while having been treated in past Cahn-Hilliard 
and Allen-Cahn approaches [8], adds additional numeri¬ 
cal challenges when also considering strain gradient con¬ 
tributions. These challenges have been overcome in this 
work as demonstrated by our three-dimensional numeri¬ 
cal examples. 

Mechano-chemical spinodal decomposition as de¬ 
scribed here can be expected in solids forming high 
temperature phases that exhibit dynamical instabili¬ 
ties at low temperature. An accumulating body of 
first-principles calculations of Born-Oppenheimer sur¬ 
faces have shown that many high temperature phases 
are dynamically unstable at low temperature [ 26 , 27 ], be¬ 
coming stable at high temperature through large anhar- 
monic vibrational excitations [ 28 - 31 , 33 ], usually through 
a second order phase transition. While such instabilities 
are frequently dominated by phonon modes describing 
an internal shuffle, a subset of chemistries become dy¬ 
namically unstable at low temperature with respect to 
phonon modes that break the symmetry of the crystal 
unit cell [ 30 , 33 ], an instability that can be described phe¬ 
nomenologically with strain order parameters. If, upon 
alloying such compounds, the high symmetry phase be¬ 
comes stable by passing through a two-phase region, its 
free energy surface will resemble that of Figures 2 b and 3 , 


and thereby make the solid susceptible to the mechano- 
chemical spinodal decomposition mechanism described 
here. 

While first-principles and experimental evidence sug¬ 
gests mechano-chemical spinodal decomposition should 
occur in ZrH2-2c5 we expect it to occur in a wide range 
of other chemistries as well. One possible example, as 
described in the introduction, is the decomposition of 
cubic yttria stabilized zirconia [ 34 , 35 ] upon quenching 
from the high temperature cubic phase into a two-phase 
region separating a low-Y tetragonal phase from a Y-rich 
cubic phase. Past treatments of this transformation [ 5 - 7 ] 
relied on non-conserved order parameters to distinguish 
the different tetragonal variants from each other and from 
the cubic parent phase. The elastic part of the free en¬ 
ergy density function was parameterized using linearized 
elasticity and infinitesimal strains, as other authors have 
also done [ 39 , 40 ]. The chemical part of the free energy 
density was assumed to have a negative curvature as a 
function of the non-conserved order parameter at Zr-rich 
compositions, but made up of convex (and quadratic) 
potentials with respect to strain. We reiterate that such 
a treatment rests on the implicit assumption that inter¬ 
nal shuffles drive the cubic to tetragonal transformation, 
while the free energy as a function of strain remains con¬ 
vex for all relevant values of the non-conserved order pa¬ 
rameters. 

The possibility also exists that a coarse-grained free 
energy density of cubic Zr02 exhibits negative curva¬ 
tures with respect to strain below the cubic to tetragonal 
transition temperature making the formalism developed 
here an accurate description of decomposition reactions 
of yttria stabilized zirconia. While the precise mecha¬ 
nism and nature of the cubic to tetragonal transition of 
pure Zr 02 remains to be resolved [ 29 ], first-principles 
calculations predict that the cubic form of Zr02 is dy¬ 
namically unstable with respect to transformation to the 
tetragonal variant [ 36 ]. A rigorous statistical mechan¬ 
ics treatment [ 28 , 33 ] is required to determine whether 
the cubic to tetragonal transition of pure Zr02 is accom¬ 
panied by a change in the sign of the curvature of the 
free energy density with respect to 62 and 63. If this 
proves to be the case, yttria stabilized zirconia should 
also be susceptible to mechano-chemical spinodal decom¬ 
position upon quenching, consistent with the coherent 
spinodal microstructures between tetragonal and cubic 
phases observed in single crystal regions of quenched cu- 
bic Zri_cYa; 02 _c /2 [ 35 ]. 

A mechano-chemical spinodal may also play a role 
in a variety of important electrode materials for Li- 
ion batteries, and intercalation compounds considered 
for two-dimensional nano-electronics. These include cu¬ 
bic LiMn204 transforming to tetragonal Li2Mn204 [ 37 ]. 
While most mechano-chemical spinodal transitions will 
occur in three dimensions, many of the qualitative fea¬ 
tures of these transitions are more conveniently illus¬ 
trated in two dimensions. Our two-dimensional studies 
should also prove relevant to understanding mechano- 
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chemical phase transformations in two-dimensional lay¬ 
ered materials for nano electronics. Materials such as 
TaS2, are susceptible to Peierls instabilities upon varia¬ 
tion of the composition of adsorbed or intercalated guest 
species that donate to or extract electrons from the sheet- 
like host [ 38 ]. Our phenomenological treatment by intro¬ 
duction of the concept of a mechano-chemical spinodal, 
coupled with gradient stabilization of the ensuing non- 
convex free energy in strain-composition space, thus of¬ 
fers a framework with potential for extension to a wide 
range of phase transformation phenomena. 

With regard to the fundamental thermodynamic un¬ 
derpinnings of analogous processes, the phenomenologi¬ 
cal model introduced here can also be used to describe 
temperature driven martensitic transformations. The 
composition, c, would then be analogous to the internal 
energy density, ^ ^ while the chemical potential, //, would 
be analogous to the temperature, T. Due to the presence 
of temperature gradients throughout the solid, though, 
the starting point must be a formulation of the entropy 
density, (as opposed to the Helmholtz free energy) as 
a functional of spatially varying ^ and displacements, u. 
In analogy with Equation ( 2 ), the total entropy can be 
written as a volume integral over a homogeneous entropy 
density that depends on the local internal energy density 
and strain ^("^,6), as well as a non-uniform entropy 
density contribution expressed in terms of gradients in 
internal energy density and strain Ve). Varia¬ 

tional maximization of the entropy will yield mechani¬ 
cal equilibrium equations as well as an expression for the 
temperature T (strictly speaking, for its reciprocal, 1 /T), 
similar to Equation ( 4 ) for the chemical potential. Due 
to the presence of gradient terms, and Ve, the tem¬ 
perature will not only be a function of the local internal 
energy density and strain, but also gradients in these field 
variables. As with the diffusion problem treated here, 
heat flow can be described with an Onsager expression 
relating the heat flux to a gradient in temperature. The 
possibility exists that the homogeneous entropy ^ ex¬ 
hibits instabilities with respect to internal energy density 
^ ^ allowing for spinodal decomposition with respect to 
the redistribution of in a manner that is similar to the 
well understood problem of spinodal decomposition with 
respect to composition. The treatment introduced here 
can therefore describe (upon replacing c with ^ and /i 
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I. THE EXISTENCE OF A FREE ENERGY 
DENSITY FUNCTION THAT IS NON-CONVEX 
WITH RESPECT TO COMPOSITION 


Consider the experimentally derived equilibrium phase 
diagram for transition metal hydrides [1] shown schemat¬ 
ically in Figure 2 c of the main text. The single cubic 
phase, stable for all composition and strains at high tem¬ 
perature, only remains stable at high composition when 
quenched to low temperature. At low composition, the 
tetragonal phase is stable at low temperature. Also re¬ 
call that first principle calculations have shown that the 
tetragonal phase is unstable with respect to strain and 
decomposes into three energetically equivalent variants 
[ 2 ]. The homogeneous free energy density therefore is 
symmetric with respect to zero strain. The two-phase 
coexistence region implies that the free energy density 
surface, ^ in strain-composition, (e, c) space admits a 
common tangent hyperplane. For clarity, let us consider 
the zero strain hyperplane in this space. Its projection 
appears in Figure SI. For a point in the solid with com¬ 
position c*, the common tangent construction leads to 
compositions Cq, and C/3, respectively of the tetragonal 
and cubic phases, with 


dc 





(1) 


Furthermore, the stability of the cubic and tetragonal 
phases with respect to composition imply that 


52^ 


dc^ 


> 0 , 


52^ 


dc^ 


> 0 . 


C/3 


( 2 ) 


Recognizing that C/3 > Cq, in Figure 2 c of the main text, 
it follows from Equation ( 2 ) that there exists c+ > Cq, 
and < Cj^ such that 


dc ^ dc 

ct Coc 


and 


d^ 

dc 
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d^ 

dc 


( 3 ) 




But then, from Equations ( 1 ) and ( 3 ), it follows that, for 
a smooth ^ that is at least twice differentiable in the 
interval (cc,C/3), there must exist at least one value c^ 
with c+ < satisfying 




d(? 


Cl 


< 0 . 


( 4 ) 


Therefore there is at least one intermediate value of com¬ 
position in the two-phase coexistence region at which the 
free energy density is non-convex with respect to compo¬ 
sition. This result is simply a consequence of the mean 
value theorem of calculus. 



EIG. SI: Figure SI: Projection of the free energy on 
to the zero strain hyperplane. 


H. FINITE STRAIN KINEMATICS 

In the interest of clarity and consistency with the main 
text, the Euclidean basis of our Cartesian coordinate sys¬ 
tem is denoted by % = 1 ,... 3 , • Xj = This 

constitutes a mild abuse of notation, because reference 
positions have also been denoted by Xi,X2,X3. The 
high symmetry, cubic reference configuration of the crys¬ 
tal is denoted by U, on which material points are labelled 
by their position vectors X = XiXi + X2X2 + X3X3. 
The deformation map from U to the current configuration 
of the crystal is = X + u, where the displace¬ 

ment field u was introduced in the main text. The de¬ 
formation gradient tensor is F = d^pjdX = 1 -\-du/dX. 
Using lower case indices to denote the components of 
vectors and tensors on the current configuration, and up¬ 
per case indices for those on the reference configuration, 
we have the coordinate expression Fij = difijdXj = 
5ij + duijdXj. 

The Green-Lagrange strain tensor is F = —1), 

where 1 is the second-order isotropic tensor. We re¬ 
peat this relation in coordinate notation from the main 
text: Eij = \{FiiFij — 5 ij). The polar decomposi¬ 
tion theorem states that the deformation gradient can 
be decomposed multiplicatively as F = Ft/, where R 
is an (orthogonal) rotation tensor belonging to the space 
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SO( 3 ). It therefore satisfies R = 1 . The other ten¬ 
sor in this decomposition, t/, is symmetric and positive- 
definite. From the polar decomposition it follows that E 
as defined above is invariant to rotation of the current 
configuration. The requirement of material frame invari¬ 
ance dictates that the dependence of the free energy den¬ 
sity on elasticity must be as a function of E. Such a free 
energy density function is invariant to rotations on the 
current configuration. 


III. THE CHEMICAL POTENTIAL: 
VARIATIONAL DERIVATION 


Since the crystal is, in general, far from chemical 
equilibrium we follow the standard treatment of non¬ 
equilibrium thermodynamics to model the chemistry. In 
order to obtain a relation for the chemical potential we 
consider variations on the composition field of the form 
Ce := c -h where ^ is a function such that, on the 
Dirichlet boundary, Tc C T, we have c = c and ^ = 0 . 

S dTl[c + e^,u\ 

( 5 ) 

Carrying out the differentiation with respect to 5 and 
integration by parts yields 

1 ^ _ 

+ E 

a ,(3 

+ E (w • - V • (0“Ve„)^ j dV 

k" Vc + Vea j • n d5, (6) 

a / 

where n is the outward normal vector to T. At chemical 
equilibrium we have 511 / 5 c = 0, yielding two conditions: 



/' 


. s ^ 9 k 

—-V • (k Vc) + Vc • Vc 

oc oc 

1 ^ _ 


a ,(3 


dd 


E(W-(^Ve«)-V.(rVe«) 


/' 


K"Vc+^6>“Veo 


n = 0 


dV = 0 

( 7 ) 

(8) 


Standard variational arguments lead to the conclusion 
that, for chemical equilibrium, the quantities within the 
brackets must vanish on the corresponding domains Q 
and T. The expression contained within brackets in the 
integrand of Equation [ 7 ] is the chemical potential /r. In 
the non-equilibrium setting that is of interest here, we 
continue to be guided by the above treatment and use 
the same definition for the chemical potential [ 3 ]: 


/i -V • (k®Vc) + Vc • vrVc 

oc oc 


+ y] ■ 


9 c 


Ve /3 


+ E ( W • - V • (r Ve„) j . 

rv / 


( 9 ) 


Variants of Equation [8] appear in a proper variational 
derivation of chemical equilibrium—a detail that is some¬ 
times overlooked in phase field treatments. We will im¬ 
pose chemical equilibrium of the boundary T, by requir¬ 
ing the term in brackets to vanish. 

With the standard, phenomenological relation J = 
—L{c^e)ViJi for the flux, the strong form of the governing 
PDE of mass balance follows: 


—J n = J 
(K"Vc+^6>“Ve„)-n = 0 

a 

c{X)=co{X) 


in Q 
on T 
on T 

at t = 0 (10) 


Note that we have assumed that there is no Dirichlet 
boundary; i.e. Tc = 0 , in stating the final strong form 
of this PDE. The conjugate boundary condition on the 
mass flux holds on the entire boundary T, in Equation 

[102] . Additionally, we have the boundary condition aris¬ 
ing from requiring chemical equilibrium on T, in Equation 

[103] . Using standard variational arguments, the weak 
form of mass balance can be obtained by starting from 
[ 10 ]. Einally, we have the initial condition on composi¬ 
tion, Equation [IO4]. The substitution of Equation [ 9 ] in 
J = —LV/r, followed by this phenomenological flux re¬ 
lation’s substitution in (lOi) shows that this is a fourth 
order PDE in space. 


A. On mass transport posed in the reference 
configuration 

We note that the mathematical formulation of diffusion 
is commonly stated in the current (deformed) configura¬ 
tion, which is the true state of the solid. However, for 
crystalline solids (with negligible extended defects) it is 
more convenient to define Onsager transport coefficients 
and evaluate composition gradients on the reference (un¬ 
deformed) configuration, here denoted by U. Diffusion in 
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the crystalline solids results from atomic hops between 
well-defined crystal sites, which can always be mapped 
onto a reference crystal structure. Increasingly, mobility 
coefficients in multi-component solids are calculated from 
first principles through the evaluation of Kubo-Green ex¬ 
pressions in kinetic Monte Carlo simulations on the refer¬ 
ence crystal in configuration Q. Since the mobility tensor 
is reported on it is most convenient to adopt the La- 
grangian description for diffusion over this configuration 
in single crystalline solids. 

Our numerical implementation is based on the weak 
form, whose derivation we begin by first reintroducing 
now in the guise of a weighting function lying in the 
function space ^ G We seek functions 

c G = {c G such that 

+ j ^ ■ {LVOV ■ (KVc)dy 

iJdS - j{LVO nV- {K.Vc)dS = 0 
r r 

(11) 



Equation [11] is obtained by multiplying ^ into [lOi], in¬ 
tegrating by parts twice, and using the Neumann bound¬ 
ary condition [IO 2 ]. Note that the higher-order Dirich- 
let boundary condition (IO 3 ) has not been built into the 
function spaces as is commonly done in weak forms. In¬ 
stead, it is imposed weakly using the method of Nitsche 
[4-6] and extending the weak form to 




'dt 


fd^\ 


y V • (ivev • {KVc)dV 


- J ^JdS - J {LVO nV- {KVc)dS 

r r 

- y V • (kVO {LVc) ■ ndS 
r 

+ J rV^ • nVc • nd5' = 0 ( 12 ) 


Equation [12] is based on the coupling tensor 0*^ = 0 
for all a, and for independent of c, corresponding 

to the numerical examples presented in the main text. 
This weak form maintains consistency. Note that the 
second-to-last term that has been added over the form 
in Equation [11] is symmetric with the term preceding 
it. This formulation shows better numerical performance 
than the also-consistent formulation in which this addi¬ 
tional term is preceded by a negative sign. Einally, r is a 


stabilization parameter that is inversely proportional to 
the element size. 


IV. MECHANICAL EQUILIBRIUM 

As explained in the main text, since elastic wave prop¬ 
agation is orders of magnitude faster than the kinetics 
of mass transport, it can be assumed that mechanical 
equilibrium is attained instantaneously for the prevail¬ 
ing composition field. To obtain the relevant equilibrium 
equations, we introduce variations on the displacement 
field: = u ^ sw, such that on the Dirichlet boundary 

Tu-, the fields satisfy Ui = Ui and Wi = 0. We construct 
the first variation of 11 [c, u] with respect to u. 



where are obtained by first replacing u by in the 
kinematic relations for E discussed above, which are 
then substituted in the relations for the reparameterized 
strains e. We recall the constitutive relations for the first 
Piola-Kirchhoff stress and the higher-order stress, respec¬ 
tively: 


_ -h dca 

^ dca dFij ^ dea,i dFij 

a a ’ 


BijK = 

a 


dcaj dFij^K 


(14) 

(15) 


Using these constitutive relations. Equation [13] can be 
manipulated to yield the Euler-Lagrange equation, which 
is also the weak form of mechanical equilibrium for a 
strain gradient elastic material, on Q in terms of P and 
B. We retain coordinate notation in the interest of trans¬ 
parency: 


J {PijWi^j + BijKWi^jx) dV - J WiFi d5' = 0 

(16) 

We introduce the normal and surface gradient operators, 
and V^, defined by 

• n, and = \/^- (V^V^)n. (17) 

Starting from Equation [16] and applying integration by 
parts several times we have the strong form of mechani¬ 
cal equilibrium for a strain gradient elastic material. See 
Rudraraju et al. [7] for detailed steps. The results of this 
variational treatment are analogous to those of Toupin 
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[ 8 ], except that it is carried out over the reference con¬ 
figuration of the crystal, rather than its current configu¬ 
ration. 


PiJ,J — BijK,JK 

= 0 

in 

n 

Ui 

= Ui 

on 


Pij'^j - '^^BijKriK'nj — 2\/j{BijK)nK 




-BiJx'^j'V^j^ + {bLLnjTlK - bjK)BijK 

= Ti 

on 


BiJKnjTlK 

= 0 

on 

r 

\n^jnKBijK\ 

= 0 

on 



Here, F = F^. U F^^ for i = 1,2,3, is the decomposi¬ 
tion of the boundary surface into a subset with Dirichlet 
boundary condition on displacement component Ui and 
Neumann boundary condition on traction component T^, 
respectively. Finally, T is a smooth boundary edge on 
which a jump in higher-order stress traction may arise. 
The fourth-order nature of the governing partial differ¬ 
ential equation emerges on evaluating [14-15] and sub¬ 
stituting in [18i]. 

The Dirichlet boundary condition in [I 82 ] has the same 
form as for conventional elasticity. However, its dual 
Neumann boundary condition, [I 83 ] is notably more com¬ 
plex than its conventional counterpart, which would have 
only the first term on the left hand-side. Equation [I 84 ] 
is the higher-order Neumann boundary condition on the 
higher-order stress traction. The physical interpretation 
of this boundary condition in its homogeneous form is 
that the boundary has no mechanism to impose a mo¬ 
ment on bonds at the atomic scale. Equation [I 85 ] spec¬ 
ifies that there is no discontinuity of higher-order stress 
traction across edges. The weak form of mechanical equi¬ 
librium has already been encountered in [16]. Eor com¬ 
pleteness we specify functional spaces: Eind u lying in 
= {u G Ui = Ui on F^^.} such that given 

functions Ti on Ft^ for i = 1,2,3 and tensors P, B de¬ 
fined by Equations [14-15], Equation [16] is satisfied for 
all w lying in y = {w G \ Wi = 0 on F^^.}. 


V. ISOGEOMETRIC ANALYSIS FOR 
HIGH-ORDER PDES 

Isogeomeric Analysis (IGA) is a mesh-based numerical 
method with NURBS (Non-Uniform Rational B-Splines) 
basis functions. The NURBS basis functions have many 
desirable properties. They are partitions of unity with 
compact support, satisfy affine covariance and provide 
certain advantages over Lagrange polynomial basis func¬ 
tions, which are the mainstay of Galerkin finite element 
methods. These advantages include the imposition of 
higher-order, C’^-continuity, positivity, convex hull prop¬ 
erties, and being variation diminishing. For a detailed 
discussion of the NURBS basis and IGA, interested 
readers are referred to Cottrell et al. [9]. However, 
we briefly present the construction of the basis functions. 


The building blocks of the NURBS are univariate B- 
spline functions that are defined as follows: Consider pos¬ 
itive integers p and n, and a non-decreasing sequence of 
values X = [^ 1 ,^ 2 , ••••,-^n+p+i], where p is the polynomial 
order, n is the number of basis functions, the ^i are co¬ 
ordinates in the parametric space referred to as knots 
(equivalent to nodes in FEM) and y is the knot vector. 
The B-spline basis functions Pi,p(0 defined starting 
with the zeroth order basis functions 


BiAO 


1 if6<^ <6+1, 

0 otherwise 


(19) 


and using the Cox-de Boor recursive formula for p > 1 

[ 10 ] 

BiAO = 

Si+p si Si+p+1 si+1 

( 20 ) 


The knot vector divides the parametric space into inter¬ 
vals referred to as knot spans (equivalent to elements in 
EEM). A B-spline basis function is C^-continuous in¬ 
side knot spans and -continuous at the knots. If an 
interior knot value repeats, it is referred to as a multi¬ 
ple knot. At a knot of multiplicity /c, the continuity is 
(jp-k ^ The boundary value problems in the main text 
consider only simple geometries, for which B-spline basis 
functions have sufficed. However, we outline the exten¬ 
sion to NURBS basis functions for the sake of complete¬ 
ness, noting that the numerical formulation as presented 
is valid for any single-patch NURBS geometry. Using 
a quadratic B-spline basis (Eigure S 2 ), a C^-continuous 
one dimensional NURBS basis can be constructed. 


v(e 


BiAOwi 

E7ABiA0n 


( 21 ) 


where Wi are the weights associated with each of the 
B-spline functions. In higher-dimensions, NURBS ba¬ 
sis functions are constructed as a tensor product of the 
one dimensional basis functions: 


II 

BiACjBjAAwij 

(2D) 

Bibl ^1)2 

E E b^aOBjAvHj 

i=l j=l 


Bi^2{0Bj,2{r])Bk,2{C)'^ijk 

nti B,h2 nhs 

- (3D) 


E E E BA2{0Bj,2{v)BkA0mjk 

i=ij=i k=i 

( 22 ) 


VI. NUMERICAL EXAMPLES 

Recall that in two dimensions, the free energy den¬ 
sity function depends only on c, 61,62 and Cq. In this 
first communication we wish to focus only on broad phe¬ 
nomenology, for which reason we choose the simplest ver¬ 
sions of the free energy density function. Apart from 
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FIG. S2: Figure S2: Quadratic (p=2) B-spline basis 
constructed from the knot vector x = [0, 0 , 0 , 1/6, 
1 / 3 , 1 / 2 , 2 / 3 , 5 / 6 , 1 , 1 , 1 ]. 


the dependence on c, which is essential for the mechano- 
chemical spinodal, the coefficients are constant. We have 
introduced the coefficients dc and de, which control the 
depths of the free energy wells and 5 e, which determines 
the strain minimum in terms of 62. Figure 2 of the paper 
depicts this construction. We also reduce k and 7*^^ to 
isotropic forms: k = and set the 

composition gradient-strain gradient coupling coefficients 
9^ = 0. While the fullest possible complexity of the cou¬ 
pling is not revealed by these simplifications, the aim 
here is to present the essential physics that is universal 
to mechano-chemical spinodal decomposition, postpon¬ 
ing details of the more complex couplings and tensorial 
forms to future communications, where they will be de¬ 
rived for specific materials systems. Accordingly, in two 
dimensions, we have: 

^ (c, ei, 62, cq) = 16 dcC^ — 32dcC^ + 16 dcC^ 

+ '" 4(4 + 4 ) 

+ ^^2 + ( 2 c - 1)^62 ( 23 ) 

^(Vc, Ve2) ' Ag/e Ve2 ( 24 ) 


where dc = 2.0, dg = 0.1, Sg = 0.1, k = 10“^ and 
Ag = 10 “^. The problem domain is a square of dimen¬ 
sions 0.01 X 0 . 01 . 


function used for the three-dimensional examples is: 

^ (c, e) =16dgc'^ — 32dgC^ + 16dgC^ 

^ + 64 -h 65 -h el) 

+ ^y2+44 

+ ( 2 c - l)^(e2 + 63) ( 25 ) 

^ (c, e, Vc, Ve) Vc • k\/c + ^ Ve2 • Ag/^Ve2 

+ • Ag/gVes ( 26 ) 

where dg = 2.0, dg = 100.0, <Sg = 0.1, k = 10“^ and 
Ag = 1 . 0 . The problem domain is a unit cube. The strain 
gradient energy coefficient tensors have been chosen to be 
isotropic, = Ag/gd"^d/j, and the three-dimensional 
formulation of the main text has been simplified here by 
setting the composition gradient-strain gradient coupling 
coefficients Ofj = 0. 

VII. LIST OF VIDEOS 



FIG. S 3 : Movie Sl[MovieSl.mp 4 ]: Evolution of 
microstructure of the two-dimensional solid with 
outward flux corresponding to the plot for = 0.1 
considered in Figure 5 a of the main text. 


In keeping with the above-explained simplifications 
aimed at presenting the essential physics of mechano- 
chemical spinodal decomposition, the free energy density 
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FIG. S4: Movie S2[MovieS2.mp4]: Evolution of 
microstructure of the two-dimensional solid with 
outward flux corresponding to the plot for = 0.01 
considered in Figure 5a of the main text. 



FIG. S6: Movie S4[MovieS4.mp4]: Evolution of 
microstructure of the two-dimensional solid due to 
quenching from an initial composition of c = 0.46. This 
corresponds to the plot for = 0.01 considered in 
Figure 5b of the main text. 
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FIG. S5: Movie S3[MovieS3.mp4]: Evolution of 
microstructure of the two-dimensional solid due to 
quenching from an initial composition of c = 0.46. This 
corresponds to the plot for /^ = 0.1 considered in Figure 
5b of the main text. 



FIG. S7: Movie S5[MovieS5.mp4]: Evolution of 
material points on the homogeneous free energy density 
surface, ^(c, ei, 62 , ce) for the two-dimensional 
formulation, restricted to the { 0 , 62 } sub-space. The 
microstructure that forms is similar to that in Movie SI. 
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FIG. S 8 : Movie S6[MovieS6.mp4]: Zoomed-in view 
with distorted mesh lines showing the evolution of an 
initially square crystal lattice (green) into rectangular 
variants (red/blue) as the chemistry-driven cubic to 
rectangular transformation occurs, followed by the 
splitting of the variants into a twinned structure. 



FIG. Sll: Movie S9 [MovieS9.mp4]: Evolution of a 
cross-hatched two-dimensional microstructure. 
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FIG. S9: Movie S7 [MovieS7.mp4]: Isovolumes of the 
three interleaved tetragonal variants and the 
three-dimensional microstructure corresponding to the 
plot for /g = 1.0 considered in Figure 6 of the main text. 
For visual clarity, only a sub-domain of the simulation 
volume is shown. 



FIG. S 12 : Movie SIO [MovieS10.mp4]: Evolution of 
plate-like structures in a three-dimensional 
EIG. SIO: Movie S8 [MovieS8.mp4]: Evolution of microstructure, 

material points on the homogeneous free energy density 
surface, ^(c, e) for the three-dimensional formulation, 
restricted to the { 0 , 62 , 63 } sub-space. The 
microstructure that forms is similar to that in Movie S7. 
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